General Relativistic Ray- Tracing Method for Estimating the 
Energy and Momentum Deposition by Neutrino Pair Annihilation 

in Collapsars 

Seiji Harikae^'^'^, Kei Kotake^'^, Tomoya Takiwaki^, and Yu-ichiro Sekiguchi^ 

^1 kkotakeSth. nao.ac.jp 



(N 



Oh; 

o 

C/2 



\o 



ABSTRACT 



Bearing in mind the application to the coUapsar models of gamma-ray bursts 
W I (GRBs), we develop a numerical scheme and code for estimating the deposition 

of energy and momentum due to the neutrino pair annihilation (i/ + z/ — )■ e~ +e"^) 
in the vicinity of accretion tori around a Kerr black hole. Our code is designed 
to solve the general relativistic neutrino transfer by a ray-tracing method. To 
solve the collisional Boltzmann equation in curved spacetime, we numerically 
^ I integrate the so-called rendering equation along the null geodesies. We employ 

the Fehlberg(4,5) adaptive integrator in the Runge-Kutta method to perform the 
^ I numerical integration accurately. For the neutrino opacity, the charged- current 

i/^ I /3-processes are taken into account, which are dominant in the vicinity of the 



accretion tori. The numerical accuracy of the developed code is certificated by 



cn I several tests, in which we show comparisons with the corresponding analytic so- 

C"^ ' lutions. In order to solve the energy dependent ray-tracing transport, we propose 

Q ' that an adaptive-mesh-refinement approach, which we take for the two radiation 

angles {6, (p) and the neutrino energy, is useful to reduce the computational cost 
significantly. Based on the hydrodynamical data in our collapsar simulation, we 
/\ ' estimate the annihilation rates in a post-processing manner. Increasing the Kerr 

c^ ■ parameter from to 1, it is found that the general relativistic effect can increase 
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the local energy deposition rate by about one order of magnitude, and the net 
energy deposition rate by several tens of percents. After the accretion disk settles 
into a stationary state (typically later than ~ 9 s from the onset of gravitational 
collapse), we point out that the neutrino- heating timescale in the vicinity of the 
polar funnel region can be shorter than the dynamical timescale. Our results 
suggest the neutrino pair annihilation has a potential importance equal to the 
conventional magnetohydrodynamic mechanism for igniting the GRB fireballs. 
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Introduction 



Gamma-ray bursts (GRBs) have long attracted the attention of astrophysicists since 
their accidental discovery in 1970s. Regarding the long-duration GRBs, there have been accu- 
mulat ing observations identifying a massive stellar collapse as their origin (e.g.. lWoosley fc Bloom 
( I2OO6I ) for a review). The duration of the lon g bursts may cor respond to the accretion of de- 
bris falling into the central black hole (BH) ( iPiro et al.lll998l ). It suggests the observational 
consequence of the BH formation likewise the supernova of neutron star formation. For their 
central engines , the so-called coUapsar has received quite some interest for more than decade 
f lWooslevl[l993l : IPaczvnskil Eggsl : iMacFadven fc WooslevI Eooof ) . 



In the collapsar scenario, the central cores with significant angular momentum collapse 
into a black hole. Neutrinos emitte d from the accret io n disk heat matte r in th e polar funnel 
region to launch the GRB outflows. iPaczynskil (jl990l ): lMeszaros fc Reed (119921 ) pioneerlingly 
proposed that the energy deposition proceeds predominantly via neutrino and antineutrino 
annihilation into electron and positron (e.g., u + i? — )■ e~ + e~^, hereafter "neutrino pair 
annihilation"). In addition, it is suggested that the strong magnetic fields in the cores of 
order of 10^^ G play also an active role both for driving the magnet o-driven jets and for 
extrac t ing a significant a i nount of energy from the centra l engin e (e.g.. iBlandford fc Znajek 



(119771 ) ; JThompson et al.l (l2004f ) ; lUzdensky fc MacFadyenI (120071 ) and see references therein) . 



However, it is still controversial whether the generation of the relativistic outflows pro- 
ceeds predominantly via magnetohydrodynamic (MHD) or ne utrino- heating processes . So 
far, m u ch attention ha s been paid to the MH D 



(2004) 



processes (e.g 



Progal (120031): 



Mizuno et al. 



Lvutikovl(l2006l):lFuiimoto et al.l (l2006[):lNagataki et al.N2q07h:lMcKinnev &: Naravan 



Komissarov fc Barkovl ( 20071 ) : lBarkov fc Komissarovl ( 20081 ) : lNagatakil ( 20091 1 : lHarikae et al 



(120071 ) 
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( 120091 ) )■ A general outcome of these extensive MHD simulations is that the magneto- 
driven shock waves can blow up massive stars along the rotational axis. Those primary 
jet-like explosions are firstly at most rn ildly relativistic due to too much baryons in the 



further out (INagatakil 



central core (e.g., iTakiwaki et al.l ( l2009l )). however could be relativistic as they propagate 



Fuiimoto et al. 



Kawagoe et al. 




20091). In such a c ollaps ar environment, explosive nucleosynthesis (e.g.. 



Hiramatsu et al. 



Nagataki et al.l (20071). neutrino and gravitational- wave signals (e.g.. 



(120051)), have been also extensively studied. 



In contrast to such blossoms in the MHD studies, there have been only a few studies 
pursuing the possibility of generating jets by the energy deposition via neutrino pair annihila- 
tion. This is mainly because the neutrino emission from the accretion disk generally becomes 
highl y asphe r ical, thus demanding us to solve a multidimensional neutrino transfer problem 



(e.g., iTubbsl ( 1l978l ): IJanka fc Hillebrandtl (Il989[ )). This is still computationally very expen- 

sive, which i s also the case for the neutrino-driven supernova simul ations (see references in 

anka et al.l (120071 )). For the first time in the collapsar simulations, iMacFadyen fc Woosley 



( 1l999l ) pointed out the importance of the energy deposition via neutrino pair annihilation. 



however the energy deposition rates to the polar funnel region were adjusted by hand to 
produce jets. To our best knowledge, the fast and coUimated neutrino-heated outfiows have 
not been realized so far in the nume rical simula t ions w ithout the ar t ificial energy injection 
to the polar funnel regions (see, e.g., lAloy et al.l ( 120001 ): IZhang et al.l ( 120031 ): iMizuta fc Aloy 
( I2OO9I ) and references therein). 



Thus far, there have been reported several methods aiming to implement the neutrino 
pair annihilation into the collapsar simulations. By estimating the fiuxes and spectra of 
the neutrino emissi o n from the accretion d isk via the so-called neutrino leakage scheme, 
Ruffert et al.l ( 1l997l ): iRuffert &: Jankal (119981 ) proposed to estimate the heating rate by sum- 
ming up the contributions of th e neutrino and a i itineu trino radiation incident from all direc- 
tions. Along this prescription, iNagataki et al.l (120071 ) have estimated the neutrino heating 
rates, and included them to the hydrodynamical simulation. For reducing the computa- 
tional time, they adde d one more assumption of the optically thinness of the accretion disk 
to the prescription by iRuffert &: Jankal (119981 ). Even with this potential overestimation of 
the heati ng rates, the neutrin o-driven outfiows were not observed in their simulations. More 
recently, iDessart et al.l (120091 ) have developed a new scheme to estimate t he energy depos i- 
tion rate using the state-of-the-art, multi-angle neutrino-transport solver (lOtt et al.ll2008al ). 
They discussed the possible formation of the neutrino-driven outfiow in the pos tmerger phase 



of bin ary neutron-star coalescence. Relying on the neutrino leakage scheme, iHarikae et al. 



(120101 ) have proposed a special relativistic ray-tracing method to estimate the annihilation 
rates. Using hydrodynamical data in their collapsar simulation, they pointed out that the 
neutrino-heated outfiow might be formed in ~ 10 seconds after the initial collapse of the 
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progenitor star. 

It should be noted that all of the above schemes neglect the general relativistic (GR) 
effects for simplicity, which have bee n reported to enhance the annihilation rates signi f icantl y 
near the accreting b l ack holes (e . g.. IJaroszvnskil (119931 Il996l ): ISalmonson fc Wilsonl ( 1l999l ): 



Asano fc Fukuvamal (120011. 120001 ) : lBirkl et al.l (120071 )). Among the GR studies, the numerical 
method of iBirkl et al.l ( 120071 ) would be one of the most sophisticated one, in which a ray- 
tracing calculation is performed to follow the neutrino trajectories in a Kerr spacetime. The 
ray-tracing method has an advantage because it can straightforwardly capture important 
GR features such as the ray bending and redshift. In their scheme, the neutrino number 
flux emitted from the accretion disk (or from the neutrino spheres) is simply assumed to 
be conserved along the geodesies. In reality, the neutrino emission, absorption, and scat- 
tering should occur along the neutrino geodesies changing its neutrino distribution function 
simultaneously. Esp ecially in the abse nce of the charged- current neutrino interactions, the 
annihilation rates in iBirkl et al.l (120071 ) could be overestimated. To improve these issues, one 
has to solve the general relativistic neutrino transport equation along each ray, which we are 
to investigate in this paper. 

In this study, we present a numerical code and scheme for calculating the deposition of 
energy and momentum via neutrino pair annihilation in a Kerr spacetime, in which we solve 
the general relativistic radiative equation along the null geodesies. The charged-current /3- 
processes are taken i nto account, which are dominant in the vicinity of the accretion tori (e.g., 
Dessart et al.l (|2009[ )). With these improvements, the newly developed code would provide 
a more realistic estimation of the annihilation rates than before. We check the numerical 
accuracy of the developed code by showing several comparison with analytic solutions, some 
of which w e newly derive in th is paper. Based on the results of our long-term collapsar 
simulation ( jHarikae et al.ll2009l ). we run our new code to estimate the annihilation rate in a 
post-processing manner and discuss their implications on the dynamics of collapsars. 

This paper is organized as follows. In Section |2l we summarize the formulation of the 
general relativistic ray-tracing method for the collisional Boltzmann equation. Section |3] is 
devoted to the numerical tests. In Section HJ we estimate the annihilation rates in a post- 
processing manner using hydro dynamical data in our collapsar simulation. We summarize 
our results and discuss their implications in Section |5l 
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Fig. 1. — Schematic picture of relations of energy /momentum among three frames (LNRF, 
BLF, RF). Upper script R denotes the variables measured in RF, while upper script L 
denotes the variables measured in LNRF. 



2. Neutrino Pair Annihilation in General Relativity 

In this section, we summarize the formalism and our strategy to estimate the neutrino- 
pair- annihilation rates based on the general relativistic radiation transfer. In section 2.1, we 
summarize the method to solve the neutrino geodesies in a Kerr spacetime for the coUisionless 
Boltzmann equation. Then in section 2.2, we move on to mention how to solve the coUisional 
Boltzmann equation along the geodesies. 

We assume that the gravitational field, which leads to ray bending and redshift, is given 
by the central Kerr BH of mass M and angular momentum parameter a = J/M (where 
J is the angular momentum of the BH, and < a/M < 1), whose metric is given in the 
Boyer-Lindquist coordinates (t, r, 9, (p) by 

ds^ = Qapdx'^dx^ , 

= -a^dt^ + -fij{dx' + I3'dt){dx^ + I3^dt), (1) 

where the lapse function a, the shift vector /3* and the non- vanishing components of the 
spatial metric 7jj are given as 



a 



"X' 



P' 



-U, Jrr = -r 



A^ 



Asin^e 



(2) 



2Mr + a^ A 



r2 + a2)2-a2Asin2^ = SA + 2Mr(r2 + a2), 



where H = r'^ + a'^ cos^ 6, A = r 

and u = 2aMr/A (e.g., iMisner et al.l ( 1l973l )). Here we use G = c = 1 unit and note that 
the Latin indices {i,j) have the domain of (r, 6, (p). For later convenience, we also define the 
dimensionless angular momentum parameter of a* = a/M. 

For later convenience, we first introduce the following three frames, the Boyer-Lindquist 
frame (BLF) which is given by the center of mass system in curved space-time, the locally 
non-rotating frame (LNRF) which is given by the tetrad frame rotating with the central BH 
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to make the dragging effects vanish (i.e., e^ = with e^ being the basis of the vierbein), and 
the rest frame of fluid (RF) which is necessary to define quantities related to radiation such 
as emissivity and absorptivity. These three frames can be connected with each other by the 
tetrad and Lorentz transformations. In the following sections, the quantities measured in the 
LNRF and RF are denoted by the superscript "L" and "R", respectively. Variables in the 
BLF are denoted without any superscripts. A schematic picture between these three frames 
are illustrated in Figure [H What we finally need is the annihilation rates measured by the 
observer in the LNRF (left-end in the figure). The neutrino emissivity and absorptivity are 
naturally defined in the rest frame of fluid (RF), in which the radiation isotropy is maintained 
(right-end in the figure). The first step is to transform the variables in the RF to the ones 
in the BLF (from the right-end to left by one step) using the tetrad transformation (the 
"relation" shown in the figure with some formulae will be derived later in this section). The 
second step is to do the ray-tracing calculation from the source to the target in the global 
BLF (indicated by "Target' " in the figure). Finally the annihilation rates in the LNRF is 
given by the tetrad transformation form the BLF to the LNRF. In the following, we explain 
these procedures more in detail. 



Th e local annihilation rate in the LNRF (e.g.. lGoodman et al. I (119871 ) : lAsano fc Fukuyama 



(120001); iBirkl et all (120071 )) is written as, 

Ql{r) = 2KGl j d'pld'p]; 

X (ej;e^)(p^ + p];),fl{p];. r)f^ip];, r) 

X [1 — sin6'i, sin6'p cos ((/9,^ — (yjp) — cos^^ cos6'p] , (3) 

where fl' is the number density of neutrinos in the phase space within the solid angle of 
dQu = smdiiddydipi, in the momentum space, p\, and e\, is the momentum and energy in 
the LNRF, respectively. Those definitions are the same for antineutrino by changing the 
notation v to v. The dimensionless parameter K is written as 

^^, _. 1 + 4 sin^ ^w + 8 sin^ ^w ,,^ 

K[Ve, Ve) = , 4 

_ _ 1 - 4 sin^ 6'w + 8 sin"^ 6'w /-.x 

A (z/^, u^) = K{ur, i^t) = 7 ■ (5) 

DTT 

Here the Fermi constant is Gp = 5.29 x 10~^^cm^ MeV~^ and the Weinberg angle is sin^ ^w = 
0.23. Since the ray-tracing calculation is conveniently done in the global BLF, we transform 
the neutrino momentum in the LNRF (p^) to Pa measured in the BLF. This transformation 
is done by the tetrad transformation as. 

Pa = uj'^aP^, (6) 
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where u^^ is the transformation matrix of the Boyer-Lindquist coordinates, 

/a 

n 



W 



a 





Ijrr 















\ 







(7) 



It is noted that the annihilation rate in the BLF is given as, 

Qm = ^;Q^, (8) 

which can be readi ly implemented in th e general relativistic hydrodynamic simulations via 
V uT^" = Q^ (e.g., IShibata et al.l (120071 )). although this is beyond the scope of this paper. 

To evaluate the annihilation rates, we have yet to determine f\--. in Equation dSD- It is 
noted that the distribution function is invariant under the tetrad transformation as 



(9) 



where /^(p) is the distribution function in the BLF. fy(,y) is d etermined by the general rela- 
tivistic Boltzmann transport equation (JMisner fc Sharplll964l ) as 



dfv{v) 

dX 

D 



P 



Pfv{y) _ (dfv(v) 



Dx^ 



d\ 



coll 



_d_ 



r^^^y 



Q7-f 



d 



(10) 



:ii) 



where {dfy(^y)/dX)^^^^ represents the collision term. In the context of photon propagation 
from the accr e tion disk, there have b een extensive studies to deterni i ne the geodesi e s (e.g . 



Cartel] JlQGsh: 



Bardeen et a. 



Ranch fc Blandford 



Cadez &: Kosti 



icl J2005h~ 



Jl994h: 



(Il972f): ICunningham fc BardeenI (119731) : ICunninghanJ ( 1975 ): 



Fanton et al. 



Cadez fc Calvanil 



(119971) 
IJ2OO5I) 



Cadez et al. (119981): ICadez et al.l (120031): 



Li et all (120051 ) : lMiiller fc Camenzindl (120041 ) 



Takahashil ( 20041 . l2005l ): iTakahashi &: Watarail ( 20071 )). Since the mass of neutrinos are neg- 



ligible compared to the relevant energy-scales to affect the dynamics of collapsars (O(MeV)), 
the neutrino geodesies can be treated as that of photon and the techniques for the photon 
transfer is also appli c able t o neutrinos. To determine the null geodesies, we basically follow 
the method in IZinkI (120081 ) which utilizes the ray-tracing method. To treat the neutrino 



transport equation with the collision term, we employ the formalism developed by iLindquist 



(119661 ). In the following, we summarize the method to determine the geodesies in a Kerr 
spacetime for the coUisionless Boltzmann equation. Then in section 2.2, we present the 
method to solve the coUisional Boltzmann equation. 



2.1. Geodesies in a Kerr Geometry 



In the Boyer-Lindquist coordinates, the Lagrangian L for d escribing the geodesies of 
massless particles in the Kerr geometry (e.g., iMisner et al.l (Il973[ )) is given as 



2C 



9al3- 



dx°' dx^ 
dX~dX 

2Mr 
1 ;— 



4aMr sin^ 6 . ■ 
-t(f) 



+W^ + 



r^ + a^ + 



S 
2a'^Mr sin^ 9 



-r^ 



A 



(12) 



where overdots denote the differentiation with respect to an affine parameter A. With three 
constants of motion. 



E = -pt, 
Lz = P4„ 
C = [Llcosec'^e-a'^E'^]cos^e + pl, 



(13) 
(14) 
(15) 



one ob tains equations governing the orbital trajectory (e.g., ICarterl (119681 ): iBardeen et al. 



p 



p' 



p 



pv- 



E — uLz 

dr 

d^ 

dO 

^^S" t d^ 



^^S" t d^ 



n 






2Mr) + 2Mr^^ 



where 



7^ 

e 

V 



A 



C - [-a^E^ + Ll cosec^ 9] cos^ 9, 



(16) 
(17) 

(18) 

(19) 



(20) 
(21) 
(22) 



The integrals of motion in Equations ( 1161) - ( 1191) can be performed either numerically or 
analytically. Although the analytic solutions, if obtained, are accurate and good for reducing 
the computational cost of the ray-tracing calculation, they may be obtained only for some 
special conditions such as the motion for = (in the (r, 9) plane). Therefore we choose to 
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perform the numerical integration, and utilize the analytic solutions to test the validity of 
the numerical integration in some test problems that will be presented in section 3. 

To capture accurately the trajectory in the vicinity of the BH, a much finer resolution 
with respect to A should be taken than for the regions far distant from the BH. Therefore 
some adaptive-m esh- re fi neme nt approach is needed for accurate and efficient numerical in- 
tegratio n. As in IZinkI (|2008l ). we ch oose to employ the scaled fourth-order Runge-Kutta 
method (JFehlberg. E.I (Il970l ). see also IPapageorgiou et all (llQSSi )). which is often referred to 
as the RKF45 method. Our choice of the Fehlberg (4,5) adaptive integrator is known to be 
very useful because it is possible to estimate a truncation error, by which the adequate step- 
sizing for the Runge-Kutta integration can be determined automatically. In this method, 
the step size in each integratioi i is controlled by co mparing the residual error dj^ to a given 
criterion 6 at every k step (see iFehlberg. E.I ( Il970l ) for more detail). Here we set 6 = 10~^, 
which provides enough accuracy in tracing the ray near the BH, as will be shown in section 
3. 



2.2. Radiative Transfer in Curved Space-time 



According to iLindquistI (11966! ) . the Boltzmann equation with the collision term for 
photons is generally expressed as. 



d\ 



n{Q- Kf). 



(23) 



Here n{x) is the proper number density of the external medium with which neutrinos interact, 
and thus measured in its own local rest frame. Q{x, p) is the emission rate per particle of the 
medium {Qe), plus a further increase due to scattering (Qs), which can be therefore written 

as 



Q{x,p) = Qc{x,e^) + Qs{x,p) 



R 3{x,e^) 



4vr(eR)2 ' 
Qs{x,p) = e^de^dn{x,p')^{x;p ^p)f{x,p'), 



(24) 
(25) 

(26) 



where j is the emissivity and ^(x; p' — )■ p) is the so-called invariant phase function, describing 
the momentum transfer due to scattering, k in Equation ( l23l) is the invariant absorption 
coefficient. dQ{x,p) is the solid angle in the momentum space of p at position x. e^ is the 
neutrino energy measured in the local proper frame that is related to the quantities in the 
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BLF as 



6^ = -P^ 



= -u^a. (27) 

The formal solution of Equation fl23|) can be given as, 



/(e,fi)= / n(A")Q(A",/)e"i'v?«(^'W^')'^^'ciA", (2^ 



which is referred to as the rendering equation of the radiation transport problem (e.g., IZink 



(120081 ) ). Note that the integration with respect to A starts from a given target point (Aq) 
where the neutrino pair annihilation occurs, propagated backward to the neutrino sources 
along the geodesies. This backward ray-tracing terminates when it hits the most outer 
boundary of our computational domain or when the optical depth for each neutrino energy 
exceeds unity indicating the surface of the neutrino spheres, both of which are represented by 
Xs in Equation ( l28l) . It is noted that we set the inner boundary of the target region to be the 
surface of the ergosphere, because we consider an ideahzed situation that the energy released 
inside the ergosphere will terminate in the BH, playing no important role to energetize a 
GRB. 



In solving the rendering equation, we neglect the scattering terms Qs (Equation (l26l)). 
which is not only difficult to be treated by the ray-tracing technique but also a major un- 
dertaking in the radiative transport problem in general. The integration in the rendering 
Equation (!28|) is done explicitly along the geodesies. In doing so, we determine each in- 
tegration step by restricting the maximum change of neutrino opacity for all the neutrino 
energy-bins to be less than 10 %. By this choice, our code can safely pass some test problems 
(see section 3). 

Neglecting the energy and momentum transfer via neutrino scattering, the neutrino 
Boltzmann equation for u^. and Ue now reads, 

^ = n[a(l - /) - /^/] = n[Q, - K*fl (29) 

where the Pauli blocking term:(l — /) is now taken into accout. It is noted that the rendering 
equation is also valid in this case by replacing k in Equation (128 p with k* = (Qg + n)- As 
for the opacity sources of neutrinos (k*), electron capture on proton and nuc lei, positron 



capture on neutron, neutrino scattering w ith nucleon and nuclei, are included ( iFuUer et al. 



19851 : iTakahashi et al.lll978l : lBruenrull985l ). Here k* is estimated as n* = S [retarget ■ o"(e^)] 



with ritargct; (^i^^) bclug the target number density of each reaction and the corresponding 
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cross section, respectively. The neutrino emission illuminated from the accretion disk mainly 
comes from the optically thick region, where the charged current /3-equilibrium should be 
nearly satisfied. Hence we estimate the neutrino emissivity as Qe = K,*f^^, where /^°[= 
l/(e'^ /"^ + 1)] is the Fermi-Dirac neutrino distribution function with a vanishing chemical 
potential. 

It is noted that an adaptive-mesh-refinement (AMR) approach that we propose in this 
paper, is an another important tool for saving the computational cost of the ray-tracing 
calculation. For example, a number of rays are required for estimating the annihilation 
rates correctly in the vicinity of the accretion disk, in which the neutrino-heated outflow is 
expected to be produced. It is therefore of primary importance to do AMR with respect to 
the angular direction of rays. Secondly, the energy bin of neutrinos is better to be treated by 
AMR, because the neutrino distribution function can be more accurately determined if the 
flner energy-bins are cast for the relevant energy scales. The actual implementation procedure 
is given as follows. Given a point x, we search the maximum intensity /(e, 6, 0) among the 
neighboring points for all the direction and for all the energy bins and call it as /max- Then 
we focus on the energy bins and angular directions, which satisfy /crit(e, 0, 0) > /C/max where 
we set /C = 0.01. Only for the domain of {e,6,(j)) satisfying the condition, we cast flner 
mesh points. In the actual implementation, we perform this selecting procedure for every 
3-dimensional space, which merits not only for saving the computational costs but also for 
maintaining the good accuracy to estimate the annihilation rates. 



3. Numerical Tests 

Before applying the newly developed code to collapsars, we shall check the accuracy of 
our code. In sections 13.11 and 13.21 we show a comparison of the neutrino trajectory between 
the numerical and analytic solution, by which we check the numerical accuracy to solve the 
coUisionless Boltzmann. In section 13. 3[ we demonstrate capability of our code to capture 
the imaging around the accreting black holes, that is the so-called BH shadow problem. In 
case of the coUisional Boltzmann equation, we perform the numerical tests to reproduce the 
radiation flelds shedding from a spherical light-bulb, which will be presented in section 3.4. 



3.1. Geodesies in the (r-9) plane 
By a straightforward, albeit tedious calculation, one can obtain the w ell-known ana - 



lytic form of the null geodesies in the {r-6) plane around a Kerr BH (e.g., ICarterl ( 1l968l ): 
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Fig. 2. — Geodesies of neutrinos in the {r-6) plane for the dimensionless Kerr parameter 
a* = (left) and a* = 0.999 (right), obtained either analytically (lines) or numerically 
(points) with the two different regulation parameters {6 = 10~^ (cross) and 10^'^ (plus)) 
for the Runge-Kutta integration (e.g., section 2.1). The central black circle (quadrant) 
represents the event horizon of the BH. It is noted that the reflection at X = or Z = is 
just for a visualization. At X = for example, the rays continue to propagate left {X < 0) 
in reality. 



Bardeen et al.l ( 119721 ): ICadez et all ( 119981 ): iLi et al.l (J2005f)). Figure |2] shows the geodesies 
near the BH in the case of a* = (left) or a* = 0.999 (right), obtained either numerically 
(points) or analytically (lines). In both cases, neutrinos are initially injected from the right 
edge of the figure with different impact parameters (for different Z in the figure). They are 
shown to be dragged by the gravity of the BH, whose surface is indicated by the black line 
in the center. It is noted that the reflection at X = or Z = is just for a visualization. 
For example at X = 0, the rays keep on propagating to the left (X < 0) in reality. For the 
numerical solutions, we vary the two different parameters {S = 10~'^,10~^), which regulate 
the numerical convergence in the adaptive integrator (see section 12. ip . We find that the 
regulation parameter oi 6 = lO"'^ is sufficient to trace the trajectory in a good agreement 
with the analytic solution, which we take in the following calculations. 



3.2. Geodesies in the (r-0) plane 

Now we move on to show the geodesies in the (r-0) plane. Since the analytical so- 
lution becomes very complicated in this case, we consider a special case that is L = aE 
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Analytical 
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r^ 




Fig. 3. — Same as Figure [2] but for the numerical solution (circle) in the (r-0) plane. r_| 
indicates the position of the outer event horizon. 



( IChandrasekharl 119831 ). In this case, the evolution equations (Equations ( !T6|) - ( TT9|) ) are 
greatly simplified as 



r = ±E, 
aE 

Combining these equations, the geodesies in the (r-0) plane becomes 



dr 



± 



a 

A' 
a 



r_L — r_ 



log — 



+ 



log — -1 

r^ — r^ \ r 



(30) 
(31) 



(32) 
(33) 



where r± is the position of the event horizon. 



Figure E] is the same as Figure [21 but for the geodesies in the (r-0) plane around an 
extremely rapidly rotating BH of a* = 0.999 (note again that a* = a/M is the dimensionless 
Kerr parameter). In the following, we call the case of a* = 0.999 as an extreme Kerr for 
simplicity. As shown, our numerical integration can reproduce the analytical solution without 
visible errors. These results support that our code can trace correctly the null geodesies in 
the Kerr geometry. 
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Fig. 4. — Comparison of the neutrino images around the accreting black holes with 
a* = 0.999, obtained either from the analytic (left) or numerical (right) integration of the 
geodesies. The neutrino flux shown is for the neutrino energy of 20MeV. 



3.3. BH Shadow for Neutrinos 

In this section, we demonstrate capabihty of our code to capture the imaging around 
the accreting black holes, which is often referred to as the BH shadow problem. 

As for the neutrino sources, we assume a thin accretion disk with a Keplerian rotation 
profile. We set the mass of the BH to be 2Mq surrounded by the accretion disk, whose inner 
and outer radius are set to be the last stable orbit of the black hole (riso) and IbGM/c^ with 
the disk thickness of vr/lO (rad), respectively. The accretion disk is set to have a uniform 
density, temperature, and electron fraction of 10^^ g cm~^, 5 x lO^^K, and 0.3, respectively. 
We focus only on the electron-type neutrino in this test problem. 

Figure H] shows one example of the neutrino images around the accreting black holes 
seen from the viewing angle of 6'vicw = 72° from the spin axis of the accretion disk. No visible 
differences are seen between the two panels, in which left and right panels are obtained either 
from the analytic or numerical integration of the geodesies. This supports the validity of our 
numerical integration of the rendering equation (Equation (128|) ). 



Figure |5] shows a variety of the images seen from various viewing angles. For example, 
when we see the accretion disk from the equatorial plane (bottom right), we can observe 
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Fig. 5. — Neutrino images of the accretion disk for different viewing angles (^vicw = 0° (top 
left), 36° (top right), from 72°, 81° to 90° (from bottom left to right). It is noted that ^vicw 
is the angle measured from the spin axis of the accretion disk. The spin parameter of central 
BH is set to be a* = 0. 



neutrinos not only from the disk of the front side, but also from the opposite side because 
of the bending of the trajectory. 

Figure |6] shows the images for different neutrino energies, while the viewing angle is kept 
fixed (6'vicw = 72°). For lower energy neutrinos (such as for 5 MeV (left panel)), the disk 
luminosity is shown to be almost north-south symmetric, while it becomes highly asymmetric 
for higher energy neutrinos (such as for 40MeV (right panel)). As the neutrino energy 
becomes lower, the position of the neutrino sphere is formed deeper inside the accretion 
disk, by which we can see the regions closer to the BH (Figure [6]). Since the angular velocity 
of the Keplerin disk is larger for the distant region from the center, the deformation of 
the images due to the Doppler effects can be more remarkably seen for the high energy 
neutrinos. It is interesting to note that in the case of the maximally rotating black hole 
(bottom panels), the BH shadow becomes asymmetric even for the low energy neutrinos 
due to the frame-dragging effects (bottom two left panels). Such features for the photon 



e = 72° 
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Fig. 6. — Images of the accretion disk for different neutrino energies {Ey) and BH spin 
parameters. The side length of each plot is 60 km. 



shadow in the vicinity of massive BHs in our Galactic center, have b een considered to give 
an important information to reveal th e mass and spin of the BHs (e.g.. lTakahashi fc Watarai 
(120071 ): iNagakura fc Takahashil (120 lOf )). Although this may not be the case for GRBs due to 
their cosmological distances, the bending of neutrinos may ha ve impac t s on th e gravitational 



radiation generated by anisotropic neutrino emission (e.g., lEpsteinl (119781 ): iKotake et al. 
(J2009al )). This can be one possible extension of this study. 



3.4. Neutrino Pair Annihilation from a Spherical Neutrino Sphere 

For the coUisional Boltzmann equation in GR, it is commonly not trivial to derive 
analytic solutions for a radiative transport problem. In the following, we derive the analytic 
solution for radiation fields, shedding from a spherical light-bulb into a uniform medium 
outside. We hope that the analytic solution may be useful to check newly developed codes 
for the radiative transport in curved space. 

In the following numerical tests, the spherical neutrino sphere with a radius of 50 km is 
assumed to have its surface temperature of T = 5 MeV on which the neutrino distribution 
function takes a Fermi-Dirac shape with vanishing chemical potential. The numerical domain 
[50 km : 300 km] is covered with Ur = 100 radial mesh points. The fiducial values of 
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Fig. 7. — Comparison of the energy deposition rate from the spherical hght-bulb test (see 
text for detail) for a given angular resolution of ray-tracing calculation {ng = 32) with or 
without the AMR treatment. Note in this test that we assume the Minkowskian geometry. 



the energy and angular bins for the ray-tracing calculation are set to be {n^L,ng,n^) = 
(16, 32, 16), which we will change to see the numerical convergence. 

To find the analytic solution in Equation (|3]), we first take the most simplest case of 
a* = OjU"^ = u^ = 0. In this case, the momenta p]^ in the LNRF can be expressed in the 
BLF as 



Pc 



^^aP^, 



(34) 



where 



w^^ = 



( i 

a 

1 



V" 



llrr 









-^ 

a 






1 

/Tee 





\ 



_\_ 



(35) 



In this way, the solid angle between the two frames can be readily shown to be the same 
{d£^ = dfi). Similarly, the volume element of the phase space and the neutrino energy in 
the local rest frame can be expressed by the variables in the BLF as follows. 



3„L 



dy 



ico'.fiPofdpodn, 



(36) 
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Fig. 8. — One example describing the casting of rays with (green) or without AMR technique 
(red). In this case, the rays are cast for estimating the annihilation rate at a given point 



(seen as a convergent point of the rays) outside the ergosphere 



ergo/ 



The concentration of 



the rays is seen (green), which is helpful to correctly estimate the heating rates with reduced 
computational cost. Note in this figure that only selected rays are chosen for illustrative 
purpose. This plot is selected from the BH shadow problems (section [3.3p for the visualization 
of the AMR. 



and 



.R 



where 






1 



(37) 



u 



1 



aZ-s-oo - gj.zi.v'^) 



3^2 ■ 



(3^ 



here we define v^ = u^ ju 



iii.fi 



Inserting these results to Equation ([3]), we obtain the following analytic forms of the 
energy and momentum deposition rate respectively as. 



Q];{r) 



2cKGlC{r)E^{r)N^{r)F{r), 
2cKGU»E^{r)N^{r)G{r). 



(39) 
(40) 
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Fig. 9. — Same as Figure [7] but for the Minkowskian case with or without rotation (circle, 
cross), the Schwarzschild case with or without rotation (square, triangle). See text for detail. 
For models with rotation, we set 1/a/1 — (v^)'^ = 2. For models with the Schwarzschild 



geometry, we put a point mass of M 



1 - {v'^y 

3Mq inside the neutrino sphere. 



Here ^^y reflects the general relativistic correction to the neutrino energy as 

-9oo{R) - 933{v'y ll-2M/R-{v^y 



-9oo[r) 



1 - 2M/r 



(41) 



where v^ = r sin 6 v^, and R is the radius of the neutrino sphere. The following two quantities 
are the energy-weighted integration of the neutrino distribution function on the neutrino 
sphere (namely fu{ry,p^)) as. 
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(42) 
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Fig. 10. — Same as Figure [7]but for the numerical solution without the gravitational redshift 
(circle), without the correction due to the tetrad transformation (triangle), and the one 
including both (cross). Agreement with the analytic solution (line) can be seen when the 
two ingredients are included (See text for detail.) 



where T^r^,) is set to be 5 MeV. Finally, geometrical factors of F{r) and G{r) are given as. 



Fir) 



G{r) 



X 



Ml — sin 9y sin 6y cos {(pi, — (fp) — cos Q^ cos 6'p] d£lud£ly 



27r2 
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,l-a;)^(l+x)(8 + 9x + 3x^), 
6 
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r) l-IMjR 



(44) 
(45) 

(46) 



To emphasize the importance of AMR for our ray-tracing calculation (e.g., section 2.2), 
we first show Figure [71 in which we compare the energy deposition rates calculated with or 
without AMR treatment. A good agreement with the analytic solution can be obtained by 
utilizing the AMR technique. We take rig = 32 with AMR to be the fiducial value in the 
following test calculations. A visualization of AMR is also given in Figure |8l 

In Figure IH we compare the analytic solutions (line) with the corresponding numerical 
solutions in the following four cases; the Minkowskian case with or without rotation (circle, 
cross), the Schwarzschild case with or without rotation (square, triangle). For models with 

2. For models with the Schwarzschild geometry, we put a 



\'i\2 



rotation, we set 1/y 1 

point mass of M = SM© inside the neutrino sphere. It is noted that analytic solutions for 

each case can be readily derived from Equation ( !39l) . In all the cases, the numerical solutions 
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Fig. 11. — Hydro dynamic configuration employed in the ray-tracing calculation. This is the 
snapshot at 9.1 s after the onset o f gravitational collap se for model JO. 8 when the accretion 
disk is in a stationary state (see iHarikae et al.l ( l2010l ) for more detail). The logarithmic 
density (in g cm~^, left-half) and temperature (in K, right-half) are shown. The white solid 
line denotes the area where the density is equal to 10^^ g cm~^, representing the surface of 
the accretion disk. The central black circle (~ 4Mq) represents the inner boundary of our 
computations. 



are shown to reproduce the corresponding analytic solutions quite well. 

Here we present the test calculations to check our implementation of the two GR factors 
in Equation (I4ip . that is the gravitational redshift {—goo{R)) and the tetrad transformation 
(~5'oo('"))- On purpose, we neglect each factor one by one, and compare it to the analytical 
solution. Figure [TO] depicts the numerical solutions including both (cross) versus without the 
gravitational redshift (circle) or without the tetrad transformation (triangle). The analytic 
solutions (lines) are shown to be reproduced only when both of them are appropriately 
included. 



4. Application to Collapsar Model 



Having checked the accuracy of our code in previous sections, we a re now in a pos i tion t o 
show an application of our code in the collapsar's environment. As in iHarikae et al.l (120 lOl ). 
we estimate the annihilation rates in a post-processing manner using the hydrodynamic 
data obtained in our long-term collapsar simulations. By comparing the ne utrino-heating 
timescale to the advection timescale of material in the polar funnel regions (see IHarikae et al. 
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( 120101 ) for detail), we discuss the possibility of generating neutrino-driven outflows there. 
Paying particular attention to the GR effects on the annihilation rates, we discuss their 
possible impacts on the coUapsar dynamics. 

As for the hydrodynamic data (such as density, electron fraction, and entropy), we take 
the ones at 9.1 s after the onset of gravitational collapse for model JO. 8 (Figure [TT]) . which 
show a clear accretion-disk and BH system with the polar funnel regions along the spin axis 
of the disk. Since this model is calculated by special relativistic hydrodynamics, we project 
those data into the ones in the LNRF for the ray-tracing calculation. The position of the 
inner boundary of the computational domain is set to be 4:Mq, which mimics the event 
horizon of the BH. We set the Kerr parameter by hand as a* = for the Schwarzschild 
geometry and a* = 0.999 for the extreme Kerr geometry. 
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4.1. Effect of General Relativity on Energy and Momentum Deposition 

To clarify the GR effects, we compare the annihilation rates in the Minkowkian (M = 
0), Schwarzschild, and extreme Kerr geometry (a* = 0.999). Figure [12] shows the energy 
deposition rate Qt (contour) (Equation[3]) and the normalized momentum transfer rate Qi/Qt 
(vector) (e.g., Equation 12]) for the Minkowskian geometry (top left/top right: with/without 
special relativistic corrections), the Schwarzschild geometry (bottom left), and the extreme 
Kerr geometry (bottom right), respectively. And Figure [13] is a difference plot, which shows 
the energy deposition rate normalized by the one in the Minkowskian geometry (top left in 
Figure [nj. 

From Figure [121 it can be seen in the Minkowskian geometry (top two panels) that the 
direction of the momentum transfer is generally radially outward, while in the Schwarzschild 
and Kerr geometry (bottom two panels), the direction especially in the vicinity of the BH 
tends to direct the center as a result of the general relativistic bending. This bending effect, 
acting to suppress the outward momentum transfer, should do harm to launch the neutrino- 
driven outflow. On the other hand, it does good to the energy deposition, because it enhances 
the head-on collision especially in the polar funnel regions. 

In fact, it can be seen that the deposition rate for the Schwarzschild and Kerr geome- 
try becomes larger than for the Minkowskian geometry (Figure [T^ . In the blueish region 
that corresponds to the polar funnel region, the energy deposition rate for the extreme Kerr 
geometry is enhanced by factors compared to the Minkowskian geometry. Interestingly it 
is mentioned that the heating rate is enhanced by about one order-of-magnitude near the 
equatorial plane in the vicinity of the BH. This is because the neutrino rays are concen- 
trated there, reflecting the conical shape of the accretion disk (triangular blackish regions 
at the sides). This concentration near the equatorial plane is found to be suppressed for the 
maximally rotating BH mainly due to the frame-dragging effect. 

To see the GR effects on the net energy deposition, we calculate the total energy depo- 
sition rate, 

■itot 



Qlu = / V^QtdV, (47) 



and the one with the outgoing momentum as I Jaroszynski] (119931 ): iBirkl et al.l (120071 )) 



Qlf = / V^Qtdv 



(48) 



->o 



where the contributions with the outgoing radial component of the momentum vector (Qr) 
are counted. As shown in Table [T] the net deposition rate and efficiency for the extreme 
Kerr geometry increase up to 16% (18% for Q°^^) compared to the Minkowskian geometry. 
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Fig. 12. — Logarithmic contour of energy deposition rates Qt [erg s~^ cm"'^] and normal- 
ized vector of momentum deposition rates Qt/Qt calculated for the Minkowskian geometry 
(top left /top right:without/with special relativistic corrections), the Schwarzschild geometry 
(bottom left), and the extreme Kerr geometry (a* = 0.999) (bottom right). The spatial 
vector is visuahzed by showing the spatial velocity vector v = Qi/Qt, which is normalized 
by the speed of hght (c = 1) being represented by the arrow (top right in each panel). The 
central black circle (~ 4:Mq) represents the inner boundary of the computational domain. 
Note that the triangular regions colored by black closely coincide the surface of the accretion 
disk (e.g.. Figure [TT]) . 
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Fig. 13. — Difference plot which shows the energy deposition rate in the Schwarzschild (left) 
and extreme Kerr geometry (right), which are normalized by the one in the Minkowskian 
geometry (see top left in Figure [T2l) . 



Our results support the previous study that GR can enh ance the heating rate, and thus 



good for the formation of the neutrino-driven outflow (e.g.. lBirkl et al.l ( 120071 )). From Table 
1, the deposition rate and efficiency are barely influenced by the spin of BH. However it 
should be note d that we have included th e spin effects only i n the radiative transport. As 
pointed out by lAsano fc Fukuyamal (1200 ll ): iBirkl et al.l ( 120071 ). the spin effects, such as on 
the structure of the spacetime (i.e., the inner-most stable circular orbit becomes smaller for 
the rapidly rotating black hole) and also on the accretion disk, should be more important 
to affect the heating rate. To clarify this point, we need a hydrodynamic data based on the 
general relativistic simulations of collapsars, which we are to investigate as a sequel of this 
study. 
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Geonietry 



Q1.P [erg s- 



<5°p*[erg s- 



efficiency 



Minkowski 
Schwarzschild 
Extreme Kerr 



6.18 X 10^° 
7.15 X 10^° 
7.08 X 10^° 



5.71 X 10^° 
6.09 X 10^° 
6.13 X 10^° 



0.510 
0.590 
0.585 



Table 1: Comparison of Q];'^*, Q^ (see Equations fl47|48p ) and efficiency for the Minkowskian, 
Schwarzschild, and extreme Kerr (a* = 0.999) geometry, which corresponds to the top left, 
bottom left and right panels in Figure [121 respectively. Efficiency is evaluated as Q^°^ / L^ 
where Ly is the total neutrino luminosity. 
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4.2. Condition for Outflow Formation 



Based on the annihilation rates in the last section, we compare the two timescales in 
this section, which are the neutrino-heating timescale and the dynamical timescale. Then 
we anticipate if the neutrino-heating outflows could or could not be produced in the polar 
funnel regions. 

To trigger the neutrino-heating explosion, the neutrino-heating timescale should be 
smaller than the advection timescale, which is characterised by the free-fall timescale in 
the polar funnel regions. This condition is akin to the condi tion of the su ccessful neutrino- 
driven explosi on in the cas e of co re-collapse supernovae (e.g.. lBethel (119901 ) and see collective 
references in iJanka et al.l (120071 )). The heating timescale is the timescale for a fluid to 
absorb the energy by the neutrino heating, comparable to the gravitational binding energy 
for making the fluid gravitationally unbound, which may be deflned as Theat = p^/Qt- Here 
$, the local gravitational potential, is t aken to be the sum of the pseudo-Newtonian potential 
and self-gravity in the flat space-time (JHarikae et al.ll2010[ ) and p is the local matter density. 



Then the dynamical timescale is deflned as Tdyn = ySTr/lGGp, where p is the average density 
at a certain radius and we take p{r) = 3M{r)/4'Kr^. 

Figure [14] depicts the ratio of the dynamical Tdyn to the heating timescales Theat for the 
Schwarzschild (left) and extreme Kerr geometry (right), showing in both cases that the ratio 
becomes greater than unity in the polar funnel regions (compare Figure [TT]) . Figure [15] shows 
the energy deposition rate along the polar axis of Figure [121 for the Minkowskian geometry 
without or with the special relativistic correction (indicated by "Newtonian" and "SR"), and 
for the Schwarzschild and extreme Kerr geometry (indicated by "GR (a* = 0)" and "GR 
(a* = 0.999)"). It is noted that the energy deposition sharply drops from the Newtonian 
to the SR case (left panel). This is the outcome of the special relativistic beaming effects. 
Since the rotational velocity of the accretion disk is perpendicular to the polar direction, the 
spec ial relativistic b eamin g effect suppresses the neutrino emission toward the polar region. 



see 



Harikae et al.l (120101 ) for more detail). When the general relativistic bending effects 
are taken into account, the deposition rate becomes larger again (see "GR (a* = 0)" and 
"GR (a* = 0.999)"). Reflecting this situation, Tdyn/Theat becomes smallest for the case with 
SR and largest for the Newtonian case (right panel of Figure [T5l) . It is important that the 



ratio in the case of the Schwarzschild and extreme Kerr geometry, which do reflect nature in 
the coUapsar's environment, becomes larger than unity inside 100 km in the vicinity of the 
rotational axis (Figure [T5l right). This indicates the possible formation of the neutrino-driven 
outflows there, if coupled to the collapsar's hydrodynamics. 
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Fig. 14. — Same as Figure [I2] but for Tdyn/Thcat (:the dynamical timescale Tdyn versus the 
heating timescale Theat) in the Schwarzschild (left) and extreme Kerr geometry (left), respec- 
tively. 
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Fig. 15. — Comparison of energy deposition rate (left) and Tdyn/Thcat along the rotational axis 
(right) between the Minkowskian geometry without or with the special relativistic correction 
(indicated by "Newtonian" and "SR"), and for the Schwarzschild and extreme Kerr geometry 
(indicated by "GR (a* = 0)" and "GR (a* = 0.999)"). 
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5. Summary and Discussion 

In the light of coUapsar models of gamma-ray bursts (GRBs), we developed a numerical 
scheme and code for estimating the deposition of energy and momentum due to the neutrino 
pair annihilation (z/ + P — )■ e~ +e~^) in the vicinity of accretion tori around a Kerr black hole. 
We designed our code to calculate the general relativistic neutrino transfer by a ray-tracing 
method. To solve the coUisional Boltzmann equation in the Kerr geometry, we numerically 
integrated the so-called rendering equation along the null geodesies. For the neutrino opacity, 
the charged-current /3-processes are taken into account, which are dominant in the vicinity of 
the accretion tori. We employed the Fehlberg(4,5) adaptive integrator in the Runge-Kutta 
method in order to perform the numerical integration accurately. We checked the numerical 
accuracy of the developed code by several tests, in which we showed comparisons with 
the corresponding analytical solutions. In order to solve the energy dependent ray-tracing 
transport, we proposed that an adaptive- mesh- refinement approach, which we took for the 
two radiation angles {6, (p) and the neutrino energy, is efficient to reduce the computational 
cost. Based on the hydrodynamical data in our collapsar simulation, we estimated the 
annihilation rates in a post-processing manner. It is found that the general relativistic effect 
can increase the local energy deposition rate by about one order of magnitude, and the 
net energy deposition rate by several tens of percents. After the accretion disk settles into 
a stationary state (typically later than ~ 9 s from the onset of gravitational collapse), we 
pointed out that the neutrino-heating timescale can be smaller than the dynamical timescale 
inside 100 km in the vicinity of the rotational axis. Our results suggest that the neutrino- 
driven outflows can possibly be launched there. 

For further investigation, we need to include several important ingredients ignored in 
this study. We plan to develop a GRMHD code for coUapsars, which is indispensable to see 
the outcome of this paper. By changing the precollapse magnetic flelds and rotation system- 
atically, we hope to clearly understand how the outflow formation in collapsars could change 
from the neutrino-driven mechanism to the MHD-driven one. The neutrino oscillation by 



the M i kheyev- Smirnov- Wolfen stein (MSW) effect (see collective references in iKotake et al. 



( 120061 ): iKawagoe et al.l (J2009[ )) could be important, albeit in much later phase than we 
considered in this paper. When the density in the polar funnel regions drops as low as 
P ^ 10^ g cm~^ later, the neutrino oscillation could operate for neutrinos traveling from the 
accretion disk to the polar funnel. If this is the case, the incoming neutrino spectra to the 
polar funnel regions and the pair annihilation rates there could be affected signiflcantly It is 
also noted that the effects of neutrino self-interaction are remained to be studie d, which has 



been attracting great attention in the theory of core-colla pse supernova e fe.g.. JDuan et al. 



( 20061)) . As in the case of core-collapse supernovae (e.g., iKotake et al.l (l2009bl ): lOtt et al. 



( l2008bl )). studies of gravitational- wave emissions from collapsars might provide us a new win- 
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dow to probe into the central engine (e.g.. lHiraniatsu et al.l ( 120051 ): ISuwa fc Murasd (120091 ) ). 



As a sequel of this work, we are planning to implement the ray-tracing calculation to the 
GRMHD simulation and clarify these issues one by one. We hope that this study takes a very 
first step towards the meeting of GR with neutrino transport, which should be indispensable 
for understanding the coUapsar engines. 
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for Computational Astrophysics, CfCA, the National Astronomical Observatory of Japan. 
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the Ministry of Education, Science and Culture of Japan (Nos. S19104006, 19540309 and 
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